home *** CD-ROM | disk | FTP | other *** search
/ Gekkan Dennou Club 147 / Gekkan Dennou Club - 2000.8 Vol. 147 (Japan).7z / Gekkan Dennou Club - 2000.8 Vol. 147 (Japan) (Track 1).bin / html / docs / complex / complex.txt < prev    next >
Encoding:
Text File  |  2000-07-11  |  8.1 KB  |  363 lines
  1. /* 複素数演算マクロ集 */
  2. /* 2000.05.21 by M.Kamada */
  3.  
  4. #include <math.h>
  5. #include <setjmp.h>
  6.  
  7. /* 演算を続行できなくなったときのためにsetjmp(cjmpenv)で環境を準備しておくこと */
  8. jmp_buf cjmpenv;
  9.  
  10. /* setjmp(cjmpenv)の返却値 */
  11. #define C_DIVIDE_BY_ZERO 1  /* 0で除算しようとした */
  12. #define C_ARG_ZERO 2  /* 0の角度を求めようとした */
  13.  
  14. #ifndef __GNUC__
  15. #error "Sorry, GCC only"
  16. #endif
  17.  
  18. #if __GNUC__==2
  19. #define complex __complex__
  20. #define Re(z) (__real__ (z))
  21. #define Im(z) (__imag__ (z))
  22. #else
  23. typedef struct {
  24.   double re;
  25.   double im;
  26. } complex;
  27. #define Re(z) ((z).re)
  28. #define Im(z) ((z).im)
  29. #endif
  30.  
  31. #define tocomplex(r,i) \
  32. ({ \
  33.   complex _z_; \
  34.   Re(_z_) = (r); \
  35.   Im(_z_) = (i); \
  36.   _z_; \
  37. })
  38.  
  39. #if __GNUC__==2
  40. #define C_1I (1.0i)
  41. #define C_2I (2.0i)
  42. #define C_3I (3.0i)
  43. #define C_4I (4.0i)
  44. #define C_5I (5.0i)
  45. #define C_6I (6.0i)
  46. #else
  47. #define C_1I tocomplex(0.0, 1.0)
  48. #define C_2I tocomplex(0.0, 2.0)
  49. #define C_3I tocomplex(0.0, 3.0)
  50. #define C_4I tocomplex(0.0, 4.0)
  51. #define C_5I tocomplex(0.0, 5.0)
  52. #define C_6I tocomplex(0.0, 6.0)
  53. #endif
  54.  
  55. #define C_2PI (6.28318530717958647692528676656)  /* 2*pi */
  56. #define C_PI (3.14159265358979323846264338328)  /* pi */
  57. #define C_PI_2 (1.57079632679489661923132169164)  /* pi/2 */
  58. #define C_PI_4 (0.78539816339744830961566084582)  /* pi/4 */
  59. #define C_1_PI (0.318309886183790671537767526745)  /* 1/pi */
  60. #define C_2_PI (0.63661977236758134307553505349)  /* 2/pi */
  61. #define C_2_SQRTPI (1.12837916709551257389615890312)  /* 2/sqrt(pi) */
  62. #define C_E (2.71828182845904523536028747135)  /* e */
  63. #define C_LN10 (2.30258509299404568401799145468)  /* ln(10) */
  64. #define C_LN2 (0.693147180559945309417232121458)  /* ln(2) */
  65. #define C_LOG10E (0.434294481903251827651128918917)  /* log10(e) */
  66. #define C_LOG2E (1.442695040888963407359924681)  /* log2(e) */
  67. #define C_SQRT1_2 (0.707106781186547524400844362105)  /* 1/sqrt(2) */
  68. #define C_SQRT2 (1.41421356237309504880168872421)  /* sqrt(2) */
  69. #define C_SQRT3 (1.73205080756887729352744634151)  /* sqrt(3) */
  70. #define C_SQRT5 (2.23606797749978969640917366873)  /* sqrt(5) */
  71. #define C_SQRT6 (2.44948974278317809819728407471)  /* sqrt(6) */
  72. #define C_SQRT7 (2.64575131106459059050161575364)  /* sqrt(7) */
  73. #define C_SQRT8 (2.82842712474619009760337744842)  /* sqrt(8) */
  74. #define C_SQRT10 (3.16227766016837933199889354443)  /* sqrt(10) */
  75.  
  76. /* xの実数部にsを加える */
  77. #if __GNUC__==2
  78. #define creadd(x,s) ((x)+(s))
  79. #else
  80. #define creadd(x,s) \
  81. ({ \
  82.   complex _z_ = (x); \
  83.   Re(_z_) += (s); \
  84.   _z_; \
  85. })
  86. #endif
  87.  
  88. /* xの虚数部にsを加える */
  89. #define cimadd(x,s) \
  90. ({ \
  91.   complex _z_ = (x); \
  92.   Im(_z_) += (s); \
  93.   _z_; \
  94. })
  95.  
  96. /* xの実数部からsを引く */
  97. #if __GNUC__==2
  98. #define cresub(x,s) ((x)-(s))
  99. #else
  100. #define cresub(x,s) \
  101. ({ \
  102.   complex _z_ = (x); \
  103.   Re(_z_) -= (s); \
  104.   _z_; \
  105. })
  106. #endif
  107.  
  108. /* xの虚数部からsを引く */
  109. #define cimsub(x,s) \
  110. ({ \
  111.   complex _z_ = (x); \
  112.   Im(_z_) -= (s); \
  113.   _z_; \
  114. })
  115.  
  116. /* xの実数s倍 */
  117. #if __GNUC__==2
  118. #define cremul(x,s) ((x)*(s))
  119. #else
  120. #define cremul(x,s) \
  121. ({ \
  122.   complex _z_ = (x); \
  123.   double _s_ = (s); \
  124.   Re(_z_) *= _s_; \
  125.   Im(_z_) *= _s_; \
  126.   _z_; \
  127. })
  128. #endif
  129.  
  130. /* xの虚数s*i倍 */
  131. #define cimmul(x,s) \
  132. ({ \
  133.   complex _z_; \
  134.   complex _x_ = (x); \
  135.   double _s_ = (s); \
  136.   Re(_z_) = Im(_x_) * (-_s_); \
  137.   Im(_z_) = Re(_x_) * _s_; \
  138.   _z_; \
  139. })
  140.  
  141. /* xのスカラーs倍 */
  142. #if __GNUC__==2
  143. #define cscal(x,s) ((x)*(s))
  144. #else
  145. #define cscal(x,s) \
  146. ({ \
  147.   complex _z_ = (x); \
  148.   double _s_ = (s); \
  149.   Re(_z_) *= _s_; \
  150.   Im(_z_) *= _s_; \
  151.   _z_; \
  152. })
  153. #endif
  154.  
  155. /* xの共役複素数 */
  156. #if __GNUC__==2
  157. #define conj(x) (~((__complex__)(x)))
  158. #else
  159. #define conj(x) \
  160. ({ \
  161.   complex _z_ = (x); \
  162.   Im(_z_) = -Im(_z_); \
  163.   _z_; \
  164. })
  165. #endif
  166.  
  167. /* xの絶対値の2乗 */
  168. #define cabs2(x) \
  169. ({ \
  170.   complex _x_ = (x); \
  171.   Re(_x_) * Re(_x_) + Im(_x_) * Im(_x_); \
  172. })
  173.  
  174. /* xの絶対値 */
  175. #define cabs(x) \
  176. ({ \
  177.   complex _x_ = (x); \
  178.   sqrt(Re(_x_) * Re(_x_) + Im(_x_) * Im(_x_)); \
  179. })
  180.  
  181. /* xの角度 */
  182. #define carg(x) \
  183. ({ \
  184.   complex _x_ = (x); \
  185.   if (Re(_x_) == 0.0 && Im(_x_) == 0.0) { \
  186.     longjmp(cjmpenv, C_ARG_ZERO); \
  187.   } \
  188.   atan2(Im(_x_), Re(_x_)); \
  189. })
  190.  
  191. /* x+y */
  192. #if __GNUC__==2
  193. #define cadd(x,y) ((x)+(y))
  194. #else
  195. #define cadd(x,y) \
  196. ({ \
  197.   complex _z_ = (x); \
  198.   complex _y_ = (y); \
  199.   Re(_z_) += Re(_y_); \
  200.   Im(_z_) += Im(_y_); \
  201.   _z_; \
  202. })
  203. #endif
  204.  
  205. /* x-y */
  206. #if __GNUC__==2
  207. #define csub(x,y) ((x)-(y))
  208. #else
  209. #define csub(x,y) \
  210. ({ \
  211.   complex _z_ = (x); \
  212.   complex _y_ = (y); \
  213.   Re(_z_) -= Re(_y_); \
  214.   Im(_z_) -= Im(_y_); \
  215.   _z_; \
  216. })
  217. #endif
  218.  
  219. /* x*y */
  220. #if __GNUC__==2
  221. #define cmul(x,y) ((x)*(y))
  222. #else
  223. #define cmul(x,y) \
  224. ({ \
  225.   complex _z_; \
  226.   complex _x_ = (x); \
  227.   complex _y_ = (y); \
  228.   Re(_z_) = Re(_x_) * Re(_y_) - Im(_x_) * Im(_y_); \
  229.   Im(_z_) = Re(_x_) * Im(_y_) + Im(_x_) * Re(_y_); \
  230.   _z_; \
  231. })
  232. #endif
  233.  
  234. /* x/y */
  235. #define cdiv(x,y) \
  236. ({ \
  237.   complex _z_; \
  238.   complex _x_ = (x); \
  239.   complex _y_ = (y); \
  240.   double _r2_; \
  241.   if ((_r2_ = cabs2(_y_)) == 0.0) { \
  242.     longjmp(cjmpenv, C_DIVIDE_BY_ZERO); \
  243.   } \
  244.   Re(_z_) = (Re(_x_) * Re(_y_) + Im(_x_) * Im(_y_)) / _r2_; \
  245.   Im(_z_) = (Im(_x_) * Re(_y_) - Re(_x_) * Im(_y_)) / _r2_; \
  246.   _z_; \
  247. })
  248.  
  249. /* xの2乗 (a+bi)^2=aa+2abi-bb */
  250. #define cpow2(x) \
  251. ({ \
  252.   complex _z_; \
  253.   complex _x_ = (x); \
  254.   Re(_z_) = Re(_x_) * Re(_x_) - Im(_x_) * Im(_x_); \
  255.   Im(_z_) = (Re(_x_) * Im(_x_)) * 2.0; \
  256.   _z_; \
  257. })
  258.  
  259. /* xの3乗 (a+bi)^3=aaa+3aabi-3abb-bbbi */
  260. #define cpow3(x) \
  261. ({ \
  262.   complex _z_; \
  263.   complex _x_ = (x); \
  264.   double _aa_, _bb_; \
  265.   _aa_ = Re(_x_) * Re(_x_); \
  266.   _bb_ = Im(_x_) * Im(_x_); \
  267.   Re(_z_) = Re(_x_) * (_aa_ - 3.0 * _bb_); /* aaa-3abb */ \
  268.   Im(_z_) = Im(_x_) * (3.0 * _aa_ - _bb_); /* 3aab-bbb */ \
  269.   _z_; \
  270. })
  271.  
  272. /* xの4乗 (a+bi)^4=aaaa+4aaabi-6aabb-4abbbi+bbbb */
  273. #define cpow4(x) \
  274. ({ \
  275.   complex _z_; \
  276.   complex _x_ = (x); \
  277.   double _aa_, _bb_, _ab_; \
  278.   _aa_ = Re(_x_) * Re(_x_); \
  279.   _bb_ = Im(_x_) * Im(_x_); \
  280.   _ab_ = Re(_x_) * Im(_x_); \
  281.   Re(_z_) = _aa_ * (_aa_ - 6.0 * _bb_) + _bb_ * _bb_; /* aaaa-6aabb+bbbb */ \
  282.   Im(_z_) = 4.0 * _ab_ * (_aa_ - _bb_); /* 4aaab-4abbb */ \
  283.   _z_; \
  284. })
  285.  
  286. /* xの5乗 (a+bi)^5=aaaaa+5aaaabi-10aaabb-10aabbbi+5abbbb+bbbbbi */
  287. #define cpow5(x) \
  288. ({ \
  289.   complex _z_; \
  290.   complex _x_ = (x); \
  291.   double _aa_, _bb_, _ab_, _aaa_, _bbb_; \
  292.   _aa_ = Re(_x_) * Re(_x_); \
  293.   _bb_ = Im(_x_) * Im(_x_); \
  294.   _ab_ = Re(_x_) * Im(_x_); \
  295.   _aaa_ = _aa_ * Re(_x_); \
  296.   _bbb_ = _bb_ * Im(_x_); \
  297.   Re(_z_) = _aaa_ * (_aa_ - 10.0 * _bb_) + 5.0 * _ab_ * _bbb_; /* aaaaa-10aaabb+5abbbb */ \
  298.   Im(_z_) = 5.0 * _aaa_ * _ab_ + _bbb_ * (_bb_ - 10.0 * _aa_); /* 5aaaab+bbbbb-10aabbb */ \
  299.   _z_; \
  300. })
  301.  
  302. /* xの6乗 (a+bi)^6=aaaaaa+6aaaaabi-15aaaabb-20aaabbbi+15aabbbb+6abbbbbi-bbbbbb */
  303. #define cpow6(x) \
  304. ({ \
  305.   complex _z_; \
  306.   complex _x_ = (x); \
  307.   double _aa_, _bb_, _ab_, _aaaa_, _bbbb_; \
  308.   _aa_ = Re(_x_) * Re(_x_); \
  309.   _bb_ = Im(_x_) * Im(_x_); \
  310.   _ab_ = Re(_x_) * Im(_x_); \
  311.   _aaaa_ = _aa_ * _aa_; \
  312.   _bbbb_ = _bb_ * _bb_; \
  313.   Re(_z_) = _aa_ * (_aa_ * (_aa_ - 15.0 * _bb_) + 15.0 * _bbbb_) - _bb_ * _bbbb_; /* aaaaaa-15aaaabb+15aabbbb-bbbbbb */ \
  314.   Im(_z_) = _ab_ * (6.0 * (_aaaa_ + _bbbb_) - 20.0 * _aa_ * _bb_); /* 6aaaaab+6abbbbb-20aaabbb */ \
  315.   _z_; \
  316. })
  317.  
  318. /* xのe乗 */
  319. #define cpow(x,e) \
  320. ({ \
  321.   complex _z_ = (x); \
  322.   double _e_ = (e); \
  323.   double _r_, _t_; \
  324.   _r_ = pow(cabs(_z_), _e_); \
  325.   _t_ = carg(_z_) * _e_; \
  326.   Re(_z_) = _r_ * cos(_t_); \
  327.   Im(_z_) = _r_ * sin(_t_); \
  328.   _z_; \
  329. })
  330.  
  331. /* xのm番目のn乗根 */
  332. #define crot(x,n,m) \
  333. ({ \
  334.   complex _z_ = (x); \
  335.   int _n_ = (n); \
  336.   double _r_, _t_; \
  337.   _r_ = pow(cabs(_z_), 1.0 / (double)_n_); \
  338.   _t_ = carg(_z_) + C_2PI * (double)(m) / (double)_n_; \
  339.   Re(_z_) = _r_ * cos(_t_); \
  340.   Im(_z_) = _r_ * sin(_t_); \
  341.   _z_; \
  342. })
  343.  
  344. /* 距離の2乗 */
  345. #define cdist2(x,y) \
  346. ({ \
  347.   complex _x_ = (x); \
  348.   complex _y_ = (y); \
  349.   double _r_ = Re(_x_) - Re(_y_); \
  350.   double _i_ = Im(_x_) - Im(_y_); \
  351.   _r_ * _r_ + _i_ * _i_; \
  352. })
  353.  
  354. /* 距離 */
  355. #define cdist(x,y) \
  356. ({ \
  357.   complex _x_ = (x); \
  358.   complex _y_ = (y); \
  359.   double _r_ = Re(_x_) - Re(_y_); \
  360.   double _i_ = Im(_x_) - Im(_y_); \
  361.   sqrt(_r_ * _r_ + _i_ * _i_); \
  362. })
  363.